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Abstract. We compare the mass functions obtained an- 
alytically, in the framework of an extended Press & 
Schechter (PS) formalism, in a previous paper to the nu- 
merical mass functions obtained in N - body simulations, 
using different algorithms to define objects in the density 
field. After discussing the properties of the algorithms, we 
show that the mass function obtained using the friend - 
of - friend algorithm reproduces best the scaling behav- 
iors predicted in the extended PS formalism. Following 
this statistical analysis, we show that it is possible in the 
framework of our extended PS formalism to reproduce the 
mass function but also, and for the first time, the initial 
statistical properties of structures and their collapse time. 
This allow to present a "coherent" picture of structure 
formation which can account for the initial, final and dy- 
namical properties of structures. 
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1. INTRODUCTION. 

This paper is the second part of a study on the influ- 
ence of non-linear dynamics on the mass function of cos- 



mic structures. In a companion paper (Audit et al. 1997, 
thereafter PI) we studied analytically the influence of the 
shear and of the tide on the mass function in a gener- 
alized Press & Schechter (1974, thereafter PS) formalism 
for a critical universe. We first showed that the mass func- 
tion can be directly related to a given dynamical model 
through a selection function which gives the probability 
that a Lagrangian fluid element with a given initial den- 
sity contrast will be in a collapsed structure at the present 
epoch. We also emphasized two invariance properties of 
the mass function in such a formalism, the well-known 
self-similarity in time, and the possibility of factorizing the 
mass function as the product of a function which depends 
only on the dynamical model and another function which 
depends only on the power spectrum of initial density fluc- 
tuations (PI). These time and spectral scaling properties 
are strong predictions of the generalized PS formalism, 



they will therefore be used to test the validity of such an 
approach. 

In the present paper, we extend our analytical work by 
a comparative analysis with numerical simulations. These 
ones have been used by many authors to test the PS ap- 
proach. The first successful attempt was made by Efs- 
tathiou et al. (1988) using however a small number of par- 
ticles. They used a rather simple and easily tractable algo- 
rithm to define collapsed objects in the numerical density 
fields, namely the Friend - of - friend algorithm (there- 
after FOF). They conclude that the standard PS mass 
function reproduces well the numerical data, although the 
number of objects with mass M ~ M* seems to be over- 
estimated (M* precedes the cut-off at large mass). Gelb 
& Bertschinger (1994) have presented a new algorithm 
(DENMAX) which brings together particles around max- 
ima of the density field. They showed that the resulting 
mass function differs from the FOF mass function. Lacey 
& Cole (1994) compared also the mass functions obtained 
with different algorithms. They used mainly FOF and an- 
other algorithm that detects spherically symmetric density 
maxima in the field (spherical overdensity, thereafter SO). 
These works emphasize that the resulting mass function 
depend strongly on the chosen algorithm to define struc- 
tures in the density field. A similar conclusion was outlined 
by Eke et al. (1996) in the case of the correlation function 
of clusters in a CDM cosmogony. 

In this paper we work then towards two goals. First, we 
study as in the previous works, the dependence on differ- 
ent definitions of structures, of the mass functions in the 
numerical simulations. In particular, we examine if the re- 
sulting mass functions satisfy the time and spectral scaling 
properties, inherent to a PS description of structure for- 
mation. Secondly, we investigate the connection between 
the actual fully non - linear dynamics of structures in 
N - body simulations and the dynamical model used to 
compute the mass function in analytical approaches. As a 
matter of fact, even if the number of objects in the field 
can be relatively well described by, for example, the stan- 
dard PS mass function, as was argued in p revious studies 
( |Efstatiou et al. 1988[ |Lacey fc Cole 199$ , this however 
does not mean that the actual dynamics of each object 
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is correctly described by the spherical model. Finally, we 
build a complete description of the mass function which 
allows both a correct statistical description of the number 
of objects generated in the numerical simulations, but also 
a right description of their initial statistical properties and 
of their dynamics. 

In the next section we present our numerical simula- 
tions of scale free power spectra and the different algo- 
rithms we used to define objects in the density field. The 
resulting mass functions and the analysis of their scaling 
properties are discussed in sect. 3. A complete description, 
which accounts for both the statistical properties and dy- 
namical origin of the mass function is exhibited in the sect. 
4. Finally we discuss our results and conclude in sect. 5. 

2. NUMERICAL SIMULATIONS 

We consider three power-law initial power spectrum 
(P(k) oc k n with n = 0,-1,-2) in a critical universe 
with a Hubble parameter h = 0.5. The simulations were 
done using a Particle - Mesh code with 128 3 particles and 
256 3 grid points to compute the gravitational force. The 
interpolation between the particules and the grid points 
was done with a Cloud-In-Cell interpolation scheme. The 
initial power spectrum were all normalized with the condi- 
tion A(i?* = 8h Mpc) — 1, where A is the r.m.s density 
fluctuation filtered by a Top-Hat window function of ra- 
dius R, and linearly extrapolated to present time. For scale 
free initial power spectra, A and the mass M are related 

by 



(i) 



where M» = (47r/3)po-R* defines a characteristic mass 
which marks the transition between linear and non - linear 
scales (A(M«) = 1). Any numerical simulations which in- 
tend to test the PS approach of structure formation should 
be able to describe correctly scales where A(M) <C 1 down 
to scales where A(M) 3> 1. Consequently, the choice of the 
simulation box size is of crucial importance. This choice is 
then a compromise between two requirements: the resolu- 
tion to describe correctly the dynamics of a single cluster 
and the power on large scales. This has lead us to choose 
a box-length of L — 80h~ 1 Mpc for n = and n = — 1 
and of L = 200ft- -1 Mpc for n = — 2 which has more power 
on large scales. This gives, with the above normalization, 
M* = 8700 particles for n = and n = -1 and M* = 560 
particles for n = —2. 

In order to consider only the structures for which the 
dynamics is properly described by the PM code, we have 
kept only structures whose virial radius (defined as the 
radius of the sphere for which the mean overdensity is 
180) are greater than 1.5 cells which is the radius at which 
the numerical force equals only half of the true force. This 
conditions gives a minimum mass of M m i n ~ 300 particles. 



This limit is very conservative but it ensures that there are 
negligible resolution effects. For each spectrum we have 
done several simulations (3 for n = 0, 4 for n = — 1 and 
5 for n — —2) to increase the statistical significance of 
our results. We obtain, for each spectrum, a few thousand 
well described (M > M min ) objects at an expansion factor 
a = 1 (present time). All the results we give in this paper 
are statistical averages over all the simulations performed 
for each power spectrum. 

The time variable was chosen in order to obtain a con- 
stant r.m.s displacement of particles, independently of the 
initial power spectrum. This can be achieve with the time 
variable p — a a where a — 2/(n + 3) (Efstathiou et al. 
1985). Each simulation was run with the same number of 
time steps (500) (using a leap-frog time integrator), which 
ensures an excellent energy conservation (see Table [1]). 

In order to study the time evolution of the mass func- 
tion we stored the results of the simulations at 6 different 
epochs. Each output corresponds to M» equal to twice the 
last output value. Table [1] gives the different simulations 
parameters we used for this study. 

2.1. Definition of objects in the density field 

To count the number of structures generated in a cos- 
mological numerical simulation, one needs to specify 
what we mean by "structure" . We have then analyzed 
the three main type of algorithms that define differ- 
ently collapsed objects in a given density field, respec- 
tively Friend - of - friend (thereafter FOF), DENMAX 



(Gelb & Bertschinger 1994), a modified version of DEN- 
MAX that we call DENEVAP, and Spherical Overdensity. 
The quality and drawbacks that one generally imputes to 
these algorithms cannot be justified by rigorous physical 
arguments. In our work we favor an algorithm, because it 
defines structures, in such a way that the resulting mass 
function satisfies the time and spectral scaling invariances, 
inherent to the PS formalism. 



2.1.1. Friend - of - Friend Algorithm 

This algorithm is one of the most used methods to define 
objects. It is based on a percolation scheme that links par- 
ticles recursively when the spatial separation between two 
particles is less than a given threshold, called the percola- 
tion length This length is expressed in units of the initial 
mean interparticular separation. This method is equiva- 
lent to selecting regions of space which are enclosed by a 
density isocontour of 1/2Z 3 . For a spherically symmetric 
object with an isothermal density profile p oc r~ 2 , this 
leads to a mean interior density p/p c = 3/2Z 3 . The choice 
I = 0.2 gives a mean density contrast of 187, which is 
the value generally assumed for the final relaxed object 
in the framework of spherical collapse. We don't however 
restrict ourselves to this value, as we will see in the fol- 
lowing. The advantages of this algorithm are its easy im- 
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#of 
runs 


(Mpc h" 1 ) 


outputs times 
(expansion factors) 


average 
AE (%) 





3 


80 


0.18 0.25 0.35 0.50 0.70 1 


0.8 


-1 


4 


80 


0.31 0.40 0.50 0.63 0.79 1 


0.9 


-2 


5 


200 


0.56 0.63 0.70 0.79 0.89 1 


0.7 



Table 1. Parameters of PM simulations used in this paper. The grid size is 256 3 , and the number of particles is 128 3 
for each run. AE is the energy conservation (in percent) for the entire simulation averaged over all the runs for a given 
spectra. 



plementation and the fact that there is only one free pa- 
rameter which can be related directly to a density isocon- 
tour. Moreover, this scheme doesn't make any assumption 
about the actual shape of objects. The main drawback 



is an overlapping proble m ( Bcrtschingcr fc Gclb 1991 



Gelb fc Bertschinger 1994 ). Indeed, small filaments can 
connect two objects which are actually different. These 
bridges of matter lead to a single massive object which 
has not collapsed, which could affect the mass function. 



2.1.2. Denmax Algorithm and Modifications 

The last remark has led Bertschinger & Gelb to develop 
a new algorithm (DENMAX) to define structures in the 
density field. DENMAX is based on a density maximum 
detection, linking together particles belonging to the same 
density peak. Particles are moved towards the peak posi- 
tion following geodesies of the density field. The surround- 
ing surface of the corresponding object is therefore not a 
density isosurface, but the col surface delineating the "in- 
fluence region" of a peak. This scheme does not suffer from 
an overlapping problem. It has however a very important 
drawback. The density maxima and the corresponding re- 
gions greatly depend on the scale at which the density field 
is computed. This means that DENMAX has a free pa- 
rameter, which is the spatial resolution of the grid used to 
compute the density field. This leads to resolution depen- 
dent results. One has to be careful when choosing the grid 
resolution. A second drawback of the method is that it de- 
fines objects without an intrinsic density criterion. Indeed, 
very low density peak are identified as objects, although 
these small density ripples are not collapsed objects. Gelb 
& Bertschinger (1994) used a additional criterion to re- 
move particles that are not gravitationally bound to the 
peak which they are associated with. Here, we rather use 
an overdensity criterion. This second algorithm, which is 
a modification to the original DENMAX, evaporates ob- 
jects which do not satisfy the following requirements: for 
each density maximum, we apply the FOF algorithm with 
a given percolation length I, and we extract the most mas- 
sive sub - halo. We then compute its center of mass, and 
its three principal axes. We then define the mean inte- 
rior density 5 of the object as the total mass divided by 
the volume of the mean quadratic ellipsoid. We vary the 
percolation length until the overdensity reaches 180. With 



this new method, low density peak are totally evaporated, 
and low density wings surrounding halos do not contribute 
to the final halo mass. For a three dimensional object with 
a complex shape , defining the mean interior density is a 
difficult task. This criterion appears to be a good com- 
promise, since we checked that in the case of a spherical 
isothermal halo, we recover the value of 180 corresponding 
to the spherical model. We refer to this modified version 
of DENMAX as DENEVAP. Note however that this new 
method still depends on the resolution used to compute 
the density field, and that the first step of DENEVAP is 
still the DENMAX algorithm. 

2.1.3. Spherical Overdensity Algorithm 

This last method is a tentative algorithm described in 
Lacey & Cole (1994) that avoids the drawbacks of the 
two previous schemes. It consists of finding spherical re- 
gions in the simulation having a certain mean overdensity 
(for example 5 — 180). We first calculate a local den- 
sity for each particle using a nearest grid point affectation 
scheme. We then sort the particles by density in decreas- 
ing order. The center of mass of a halo is first identified 
with the first particle of the list. We grow a sphere cen- 
tered on this point until the mean overdensity reaches the 
chosen value. We then calculate the new center of mass of 
the particles contained in this sphere, and we iterate until 
the distance between two consecutive center of mass is less 
than 10% of the mean interparticular distance. We finally 
remove all particles of the new halo from the list. Note 
that the halos do not depend on the method we use to 
compute the density, and to find the density maxima, as 
long as we use a high resolution grid (we use a 1024 3 grid 
in this paper). This algorithm is thus resolution indepen- 
dent. There is no overlapping problem, and there is only 
one free parameter, namely the chosen overdensity thresh- 
old. However, this scheme has the drawback of imposing 
a spherical shape on the halos. 



3. Multiplicity Function 

In this section we present the mass function obtained at 
different times, for the different power spectra and for the 
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Fig. 1. Multiplicity functions obtained for DENEVAP, SO60 and FOF-0.3 and for the different initial power spectra. 
Each curve corresponds to a different output time rescaled to a = 1. The two solid lines represent the PS mass function 
with 5 C = 1.4 (largest high-mass cut-off) and 1.686. 



algorithms presented above. When looking at these mass 
functions we have two goals. First we want to look at their 
dependence on the different methods to define structures. 
Then, we want to see if they satisfy the time self - sim- 
ilarity and the spectrum invariance that are very generic 
features of the mass function obtained through a PS - like 
approach (PI, see Appendix). 



Following Efstathiou et al. (1988), we use the multi- 
plicity function defined as the fraction of mass embedded 
in halos with a number of particles between 2 m and 2 m+1 , 
where m is the multiplicity. In these terms, the multiplic- 
ity of the characteristic non - linear mass M* between two 
consecutive outputs increases by one. In Figs. |l| and || we 
present the multiplicity functions obtained for each spec- 
trum and for each algorithm. The multiplicity functions 
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Fig. 2. Same as figure (0) but for DENMAX SO180 and FOF-0.2. 



obtained at each time are rescaled to a = 1, in order to 
test the time self - similarity of the mass function. Note 
also that the mass (in number of particles) of halos ob- 
tained in the n = — 2 simulations were multiplied by a 
factor 2.5 3 to account for the different box size. The error 
bars in both figures are computed using a Poisson noise 
estimation in each bin. We now discuss the multiplicity 
functions obtained with each algorithm. 



3.1. FOF 

As we have seen this algorithm has one free parameter, 
the percolation length I. We have first tried the value 
I = 0.2, adopted by several authors. Our results (Fig. 2) 
are in good agreement with Efstathiou et al. (1988). The 
time self similarity is very well recovered and the high 
mass cut off is well reproduced by the standard PS mul- 
tiplicity function (<5 C = 1.686, see appendix). However, as 
mentioned by Efstathiou et al., there is a mass deficit of 
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about 50% around the peak of the multiplicity function. 
This deficit suggests that the density criteria is too strin- 
gent. Therefore we have tried FOF with I = 0.3. This 
roughly corresponds to halos with an overdensity of 60. In 
this case the cut-off corresponds to a PS multiplicity func- 
tion with S c = 1.4 and the discrepancy around the peak 
is considerably reduced (Fig. 1). Since for both the per- 
colation lengths used the time self-similarity is very well 
recovered it is possible to average the multiplicity func- 
tion obtained at different times (rescaled to a = 1) to get 
a unique multiplicity function for each spectra. One can 
then check that all the multiplicity functions obtained for 
the different spectra have the same universal multiplic- 
ity function (see appendix or PI for definition). We have 
plotted on Fig. 4-A (resp.4-D) the universal multiplicity 
function for I = 0.3 (resp. I = 0.2). The "spectral scaling" 
predicted in the framework of a PS approach is extremely 
well recovered. 

For the FOF algorithm, and for both percolation 
lengths, the generic features (time and spectral scaling) 
predicted by the PS formalism are well reproduced. For 
a percolation length of 0.2 the cut-off of the multiplicity 
function is well reproduced by the spherical PS multiplic- 
ity function but in that case there is a large excess of 
"analytical" structures around M*. For a larger percola- 
tion length (I = 0.3) the entire multiplicity function is 
quite well reproduced by a PS function with 5 C = 1.4 
(even though there is still a small excess of "analyti- 
cal" structures around Af„). We have shown in PI that 
a PS-like function with S c = 1.4 results from the PS for- 
malism with a non-spherical dynamics. We will see in the 
next section if this agreement can really be explained by 
this dynamical reason. 

3.2. Denmax 

For this algorithm the free parameter is the resolution of 
the grid on which the density field is computed. Since we 
use a PM code we have chosen the same grid as the one 
used during the simulation. Using a finer grid would not 
make much sense since density maxima within a cell of 
the PM do not have any meaning. On Fig. 2 we present 
the multiplicity function obtained with DENMAX and 
rescaled to a = 1 as before. For n — the time scaling 
property is well recovered even though one starts to see 
a little shift in the small mass region and the agreement 
with a PS function with S c — 1.4 is striking. However for 
n = — 1 and n = — 2 there is a clear shift between the 
multiplicity function at different times. The time scaling 
property is not satisfied at all by the structures recog- 
nize by DENMAX. As can be seen on Fig. 4-B where we 
have plotted the universal multiplicity function of DEN- 
MAX for the three different power spectrum the spectral 
invariance is not satisfied either (We have done an average 
of the multiplicity function at different times in order to 
have a single curve even though it is not legitimate in this 



case). This fact can also be seen on Fig. 2 by looking at 
the relative position of the multiplicity function compared 
to the PS multiplicity function. The large modification 
of the DENMAX multiplicity function when one changes 
the spectrum of initial fluctuation can be understood by 
looking more closely at the structures found. In Fig. || we 
have plotted the average density versus the mass of all the 
structures found for the different spectra. It is now quite 
clear that the structures found by DENMAX greatly de- 
pend on the spectra. For n = the density wells are quite 
steep and therefore the DENMAX structures are dense. 
On the contrary for n = —2 the density wells are rather 
flat which leads to less dense structures. More generally, 
DENMAX detects as structures regions which are not nec- 
essarily dense. In order to correct this drawback we have 
developed the algorithm DENEVAP presented above. In 
Fig. |l| the multiplicity function obtained with DENEVAP 
is presented. By taking away the low density region we 
have improved the time scaling at high mass, but it is still 
very bad (worse) for intermediate and low mass. It is of 
course possible to truncate the DENMAX structures in 
many other ways but the results should not be very differ- 
ent. Beyond the problem of the structures' densities, the 
DENMAX algorithm is not very well suited to the detec- 
tion of cosmic structures because, with the grid used to 
compute the density field, it introduces a characteristic 
scale, L. This length, fixed once and for all, correspond to 
a density fluctuation, A, which is different depending on 
the time and on the spectrum considered. For this reason, 
the two scaling behaviors are broken. To recover them, 
the length L could be determined as a function of the 
time and the power spectrum in order to correspond to a 
constant A. But this artificial introduction of the scaling 
laws in the algorithm used to detect structures does not 
seem satisfactory to us. 

3.3. Spherical Overdensity 

The only free parameter of this algorithm is the den- 
sity threshold. To compare with FOF we have tested two 
thresholds: 180 (SO180) and 60 (SO60) which roughly 
correspond to FOF with a percolation length of 0.2 and 
0.3 respectively. The resulting multiplicity functions are 
presented on Figs. [I] and ||. The time self similarity is 
roughly reproduced for SO60 but very badly for SO180. 
We have therefore also looked at the spectral invariance 
for SO60. The results are much better than those of DEN- 
MAX but there is a slight shift. FOF reproduces this in- 
variance much better. For n = —2 structures are rather 
isolated from one another and therefore imposing a spher- 
ical shape is not a strong constraint. However for n = — 1 
and n = structures are more strongly correlated and 
imposing a spherical shape on an heavy cluster damages 
close-by structures. SO180 works poorly because imposing 
a spherical shape to structures with such a high density 
criteria results in the truncation of some dense parts of 
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Fig. 3. Average density of structures found by DENMAX plotted against the logarithm of their masses for the three 
different initial power spectrum. 



elongated structures. For SO60 this truncation occurs less 
and concerns region of inferior density contrast. 



4. Link with the initial conditions 

In the previous sections we have analyzed the agreement 
between the multiplicity function of numerical structures 
recognized with a given algorithm and the ones obtained 
in a PS formalism using a given dynamics. This analysis 
was only statistical and its sole concern was the number 
of structures. However, this point of view is not enough 
to have a coherent description of structure formation. If a 
numerical and an analytical multiplicity function are iden- 
tical it is then also necessary that the numerical structures 
follow the dynamics used in the analytical picture. Oth- 
erwise the agreement of the multiplicity functions is a fit 
rather than a dynamical description of structure forma- 
tion. In the PS formalism the dynamics intervene only 
through the collapse time. We have therefore decided to 
compare the collapse time of the structure predicted by 
dynamical models and the one obtained in the simulation. 
In order to give a satisfactory description of structure for- 
mation, the PS formalism should agree with the numerical 
simulations both in the initial and final conditions and also 
in the dynamical link between these two states. We will 
test the agreement in the initial conditions by comparing 
the initial distribution of shear and density contrast of 
the selected halos, for the final conditions we will use the 
multiplicity function and the dynamical coherence will be 
tested by comparing the collapse time. 

We have done this global analysis with the halos de- 
tected by the friend-of-friend because they satisfy better 
the two scaling criteria. 



Computation of the initial parameters of halos 

The first step of this analysis is to determine the initial 
properties of the halos detected in the final state of the 
numerical simulation. Once a group of particles has been 
identified as a halo of the final density field, we determine 
its initial lagrangian volume (i.e. the volume occupied by 
theses particles in the initial conditions) . We then compute 
the density field in the initial condition with a cloud in 
cell scheme using all the particles of the simulation. The 
shear which is the traceless symmetric part of the tensor 
(dvi/dxj), is then deduced from the density field using the 
relations dvi/dxj = d 2 S/dxidxj valid in the linear regime. 
The initial density and shear of a halo are then defined 
as the average of the previous fields over its lagrangian 
volume. The resulting shear tensor is then diagonalized to 
obtain the initial shear eigenvalues. 

With the previous procedure, the density contrast that 
we obtain corresponds to the value of the initial density 
field smoothed with a top-hat filter having the shape of 
the initial lagrangian volume. However, the variable used 
in the PS formalism is the initial density field smoothed 
by a spherical top-hat filter. It is therefore necessary to 
check that these two density fields have identical statistical 
properties. This has been done by comparing directly the 
probability distribution function (PDF) of both fields (for 
several halos shapes) . The maximum discrepancy between 
the two PDF's was found to be in all cases less than 5%. If 
these two smoothing methods are statistically equivalent, 
they however differ drastically on a halo-by-halo basis. 
The initial conditions attributed to a given structure will 
be very different depending on whether they are computed 
in a sphere or in the initial lagrangian volume. From an 
analytical point of view, the choice of the filter has no 
influence on the shape of the mass function for scale-free 



initial conditions (Lacey & Cole f993) 
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Fig. 4. Universal multiplicity functions obtained for FOF-0.3, DENMAX, SO60 and FOF-0.2. Each curve corresponds 
to a different power spectrum. As before, the solid lines correspond to PS mass function with 5 C = 1.4 and 1.686. 



4-2. Toward a coherent picture ? 

Now that we have both the initial properties of the struc- 
tures detected by FOF-0.2 and FOF-0.3 and their multi- 
plicity functions, we are in a position to test the dynamical 
meaning of the PS formalism as a tool to study structure 
formation. 

Analytically, structures can be identified, as we have 
stressed in PI, with regions which have collapsed either 
along their first, second or third principal axis. Before 
looking for a precise dynamical model, it is therefore nec- 
essary to determine which of these three type of collapses 
corresponds to the structures found in numerical simu- 
lations. As we have shown in PI, the derivative of the 
selection function makes a direct link between the multi- 
plicity function and the underlying dynamical model. For 
this reason, it has very specific behaviors corresponding to 
each of these three type of dynamics. Therefore, we will 
find if the structures found in the numerical simulations 



result from a first, second or third axis collapse by looking 
at the derivative of their selection function. 



From formula A. 3 (in the appendix), it is possible to 
deduce that the function P A {6) oc e'^ / 2A " dS(6, A)/dA 
gives the probability that a halo of mass A has an initial 
density contrast 5. Choosing a mass range between A and 
A + dA, it is possible to compute this initial probability 
distribution function using the results of our numerical 
simulations and then to infer the derivative of the selec- 
tion function. We plot on Fig || the derivative of the selec- 
tion function for three different values of A. The observed 
behavior (the position and the width of the peak are in- 
creasing functions of A) is characteristic of a third axis 
collapse (PI). The collapse along the third axis is slowed 
down by the effect of the shear. Therefore, on small scales 
(high A) where the shear is statistically greater, structures 
need on average a higher density contrast to collapse. We 
therefore deduce qualitatively that the underlying dynam- 
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Fig. 5. Derivative of the selection function obtained from numerical simulation data for FOF-0.3 (upper-left panel), 
FOF-0.2 (upper-right panel). The two lower panels present the same function corresponding to the analytical model 
(see text for the definition of the parameters) using the same bin size. The solid, dotted and thick lines correspond 
respectively to A =1.5, 2.5 and 3.5. 



ical model, in the framework of the PS formalism, is better 
described by the collapse of the third axis. 

Since the collapse of the first axis is a singularity, it is 
not possible follow the dynamics of the structure until the 
collapse of the third axis. However, extrapolating results 
obtained for the collapse time of the first axis, we have 
proposed in P 1 the following dynamical prescription which 
allows to compute the collapse epoch of the third axis, as 
a function of the initial density contrast S and the largest 
(and therefore positive) shear eigenvalue a 

1 - S i 1 ° V ( 2 ) 

where 5 C is a density threshold, a c is the maximum allowed 
initial shear-to-density ratio beyond which shear effects 
prevent any halo from collapsing. The last parameter, e, 
determines the global shape of the collapse epoch. From 
asymptotic arguments, we stressed in PI that e has to be 
less than unity. 

In order to recover a coherent description of structure 
formation within the PS approach, we have verified that 



a single set of parameters is able to reproduce simultane- 
ously the multiplicity function (final state), the mean and 
the variance of the function dS/dA (initial state), and 
finally the collapse time of halos. 

For FOF-0.3 (resp. FOF-0.2) we find that a good agree- 
ment is obtained with S c = 1.35, a c — 1.2 and e = 0.3 
(resp. S c — 1.4, er c = 0.6 and e = 0.7). The correspond- 
ing collapse epoch are plotted in Fig. The two collapse 
times corresponding to each algorithm are quite different. 
But they however result from the same underlying non- 
linear gravitational dynamics. The difference comes from 
the fact that the two algorithms define different regions 
in the final density field and that the initial physical con- 
ditions necessary to the formation of such region are thus 
different. The collapse time corresponding to FOF-0.2 is 
always larger. This is natural since it takes more time for 
a given halo to reach a higher density contrast. The crit- 
ical shear value is greater for FOF-0.3 than for FOF-0.2. 
This was also to be expected since, in the context of our 
dynamical model, the shear always slows down the col- 
lapse. The largest peaks of the final density field therefore 
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Fig. 7. The two upper panels show the numerical mass function obtained for FOF-0.2 (left) and FOF-0.3 (right). The 
solid curves are the mass function obtained with the analytical model presented in the text. The middle panels show 
the average (square) and dispersion (triangle) of the function dS/dA with the curve corresponding to the model. The 
lower panels shows the average and the dispersion of the collapse time found for the numerical structures with the 
analytical model. 



come from region with a low shear to density ratio. This 
result is consistent with those obtained by other authors 



(Van de Weygaert fc Babul 1994, Bernardeau 1994), but 



we analyze it from a very different point of view. The fact 
that the dense regions in the final density field have less 
initial shear does not mean that the dynamics is quasi- 



spherical, but, on the contrary, that the dynamical influ- 
ence of the shear is such that dense structure can form only 
in region where it is sufficiently low. The dynamics of the 
densest structures is a posteriori spherical (because the 
spherical collapse is the fastest) but one cannot suppose 
a priori that the dynamics is spherical everywhere. 
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0.8 



FOF-0.2 




a/6 



Fig. 6. Inverse of the collapse time (formula (g) for 8 — 1 
as a function of the shear. The solid (resp. dashed) line rep- 
resents the model corresponding to FOF-0.3 (resp. FOF- 
0.2). 



We plot in Fig. [jjthe best agreement obtained for FOF- 
0.2 and FOF-0.3 with the collapse time given by equation 
(^J) and the above parameters. We reproduce in both cases 
very accurately the multiplicity function, with a slight un- 
derestimation at the very low mass end. Note that the er- 
ror bars are only Poisson noise estimators which do not 
take into account the systematic errors coming from the 
finite resolution of the PM code such as, for example, the 
loss of dynamical accuracy in the halos of lower mass. 
Moreover, we also reproduce the mass dependency of the 
mean (< 8 >) and of the variance (E) of the initial den- 
sity contrast distribution of the halos (Fig. fj]). We have 
also plotted in Fig. [| the entire derivative of the selection 
function. 

Finally, the mean collapse epoch is recovered (a c — 
1), although the dispersion remains large for small mass 
objects. Note that for small mass halos, in addition to 
lower dynamical accuracy, there is also a systematic error 
in the computation of the initial parameters. In fact, due 
to the discreteness of the computational grid, the error 
done when averaging over the Lagrangian volume then 
increases because boundary effects become important. 

These results show altogether that our formalism is 
able to reproduce the mass function both from a statis- 
tical and from a dynamical point of view. However, each 
halo does not individually follow exactly our simple dy- 
namical prescription but in each mass range the collapse 
time of the halos is on average given by our model. For 



the sake of comparison, we plot in Fig. || the same results 
obtained with the standard PS approach, using 8 C as the 
only free parameter. Wc plot in each panel the results ob- 
tained using two different values for S c , showing that it 
is impossible to reproduce simultaneously the final and 
initial states. In other words, the standard PS approach 
gives a rough statistical indicator of the number of halos, 
but does not allow for a coherent description of structure 
formation. 



5. Conclusion 

In PI we have extended the PS formalism in order to in- 
clude the effect of the shear and the tide. We have shown 
than within the PS formalism it was possible to obtain 
very different behaviors for the mass function which are 
directly related to the underlying dynamical model. Even 
though the mass functions are quantitatively very differ- 
ent, they all satisfy, within the PS formalism, two scaling 
properties (i.e. the time self similarity and the power spec- 
trum scaling (formula (A.l)). The aim of this paper was 



to determine if the PS formalism could explain coherently 
structure formation both from a statistical and a dynam- 
ical point of view. 

In the first part of this paper we have shown, using 
numerical simulation, that it was possible to define nu- 
merical structures whose mass function obeys both the 
PS scaling laws. We have tried different algorithms which 
all lead to different mass functions. The only one which 
reproduces all the scaling laws is friend-of-friend. How- 
ever, the resulting universal mass function is still strongly 
dependent on the percolation length. This dependence of 
the mass function on the chosen algorithm shows that the 
choice of the dynamics used in the PS formalism cannot 
be unique and has to be related to the actual definition 
of the structures in the density field. In order to have 
a coherent description of structure formation, the cho- 
sen dynamical model should not only reproduce the mass 
function, but should also give a correct description of the 
regions selected in the initial density field and of their 
collapse times. We have shown, using a dynamical pre- 
scription proposed in PI, that it is possible to built such a 
coherent picture. Within this picture, the structures result 
from a collapse along their third principal axis. The col- 
lapse epoch is then a function of the density contrast and 
of the largest shear eigen- value given by formula (||) which 
has three free dynamical parameters. It is then possible to 
find a set of parameters which reproduces simultaneously 
the mass function, the initial selection function and gives 
on average the correct collapse epoch. Such a coherent 
description is totally impossible with the classical PS for- 
malism even when adjusting the density threshold. This 
inadequacy comes from the fact that the shear plays an 
important role in the dynamics of forming halos. This dy- 



12 



Non-linear-Dynamics and the mass function 



0.5 
0.4 
0.3 

0.2 
0.1 






2 

1.5 



H 1 I I I I 




H 1 I I I 



0.5 



_1 I I I I L_l_ 



_1 I I I I I I L 



_1 I I I I L_l_ 



_1 I I I I I I I 



10 



10 



Fig. 8. Same as figure ((?]) but for the classical PS formalism. The solid lines correspond to S c = 1.686. The dashed 
lines correspond to S c = 2. (left column) and S c — 1A (right column). 



namical effect explains the 50% discrepancy of structures 
found by FOF-0.2 compared to the classical PS formalism. 
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Appendix A: Appendix 

The mass function, <&(M) which gives the number den- 
sity of object of mass M, and the multiplicity function 

/u(M), which gives the fraction of mass embedded in ob- 
jects of masses between M and M+dM are simply related 
by: $(M) = ^p{M) where p is the mean density of the 
universe. In the framework of a PS approach the multi- 
plicity function can be expressed as: 



where A is the r.m.s density contrast filtered by a Top 
- Hat window function of radius R (This corresponds to 
objects of mass M = {An / '3) poR 3 ) ■ The first term depends 
only on the power spectrum of the initial density fluctua- 
tion, while the second term which we call the universal 
multiplicity function depends only on the dynamical 
prescription describing the formation of structure. Assum- 
ing spherical dynamics p takes the form: 



with S c = 1.686. More generally, the universal multiplicity 
function can be expressed as (PI) 



where C is a normalization constant. The derivative of 
the selection function, S, depends only on the dynamical 
model chosen to compute the collapse time of structures. 
All the dynamical information needed in the PS formalism 
is contained in this function. 
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